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Thermodynamic properties of the SU(n) Heisenberg model in one dimension is studied by 
means of high-temperature expansion for arbitrary n. The specific heat up to 0[(/3J) 23 ] and the 
correlation function up to 0[(/3J) 18 ] are derived with /3J being the antiferromagnetic exchange 
in units of temperature. It is found for n > 2 that the specific heat shows a shoulder in the high- 
temperature side of a peak. The origin of this structure is clarified by deriving the temperature 
dependence of the correlation function. With decreasing temperature, the short-range correlation 
with two-site periodicity develops first, and then another correlation with n-site periodicity at 
lower temperature. This behavior is in contrast to that of the inverse square interaction model, 
where the specific heat shows a single peak according to the exact solution. Our algorithm has 
an advantage that neither computational time nor memory depends on the multiplicity n per 
site; the series coefficients are obtained as explicit functions of n. 

KEYWORDS: high temperature series expansion, orbital degeneracy, Heisenberg model, SU(n) symmetry, multi- 
polar interaction, fluctuation, exchange interaction model 
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Recently, orbitally degenerate systems attract much 
attention as strongly-correlated electron systems. The 
orbital degeneracy is regarded as an internal degrees of 
freedom in addition to the spin degeneracy. The total 
number n of internal degrees of freedom thus varies de- 
pending on the system. If the number n 2 — 1 of compet- 
ing fluctuations is large, the transition temperature of 
the orbital ordering is significantly suppressed from the 
the mean- field value even in three dimensions. 1 ^ Thus 
reliable theoretical scheme beyond the mean-field theory 
should be necessary to investigate orbitally degenerate 
systems. 

As the simplest model to describe the orbitally degen- 
erate systems, the SU(n) Heisenberg model has been in- 
vestigated. In one dimension, it can be solved exactly by 
the Bethe ansatz. 2 ) However, some interesting quantities 
like the correlation function are difficult to be derived by 
the Bethe ansatz. Some critical exponents have been cal- 
culated with use of the non-abelian bosonization 3 -* and 
the conformal field theory. 4 ) Numerical methods have 
also been used to study the system. For n = 4, the spin- 
spin correlation function has been calculated using the 
density matrix renormalization group 5 ) and the quan- 
tum Monte Carlo method. 6 ) The latter has also derived 
temperature dependence of the correlation lengths and 
the entropy. Although these numerical approaches are 
powerful, they are not always convenient to investigate 
systems with large n, or in higher dimensions. 

As an alternative approach which is applicable to any 
dimension, the high-temperature expansion has provided 
significant information on variety of models. 7 ) We have 
recently developed an efficient algorithm to carry out 
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the high-temperature expansion for systems in which the 
Hamiltonian is written in terms of exchange operators. 
With use of this new algorithm, we report in this paper 
a peculiar temperature dependence of the spatial corre- 
lation, which becomes conspicuous as the number n of 
internal degrees of freedom increases. Namely, with de- 
creasing temperature, the short-range correlation with 
two-site periodicity develops first, and then another cor- 
relation with n-site periodicity at lower temperature. 
This leads to a corresponding structure in the specific 
heat. 

In our algorithm, neither computational time nor 
memory depends on the multiplicity n per site; n is just 
a parameter — the series coefficients are obtained as ex- 
plicit functions of n. Details of the algorithm as well 
as application to higher dimensional models will be re- 
ported separately. 

We consider a one-dimensional lattice of N sites with 
a periodic boundary condition. Let each site take one of 
the n colors, and denote them as |a) with a = 1, 2, • • • , n. 
We define X a P := |a)(/?| and an exchange (or permuta- 
tion) operator 



P — 



<*=! /3=1 



(1) 



Colors of sites i and j are exchanged when Pi j is applied. 
The Hamiltonian is given by 



JV-l 



N • 



(2) 



Then, the Hamiltonian has the SU(n) symmetry. We 
investigate the antiferro- interaction model with J > 0. 
When n = 2, this model is equivalent to the ordinary 



1 



2 



Letters 



Heisenberg model. The SU(4) Heisenberg model is re- 
lated to orbital- and spin-degenerate systems. The local 
four states are denoted as |+)|+], |+)|— ], |— )|+], |— )|— ], 
where |±) and |±] stand for a orbital state and a spin 
state respectively. The pseudo-spin operators are de- 
fined as t z \±) = ±±|±), t±| T ) = |±), s z \±] = ±±|±], 
s ± \^f] = |±], and then the exchange operator for n = 4 
is rewritten as 

Pi^j — ti ' tj S{ • Sj -\- 4:(ti • tj^^Si • Sj^ -\- — , (3) 

or, 

APiJ = Ti ■ Tj + <J l ■ (Tj + (Ti ■ Tj)(<Ji ■ <Jj) + 1 (4) 

using the Pauli matrices t — 2t and er = 2s. 

The high temperature expansion is performed by ex- 
panding the Boltzmann factor e~^ w in 0. To investi- 
gate the specific heat, we calculate (Pis+i) because the 
internal energy is obtained by (Ti) = NJ(Pi t i + \). Wc 
also investigate correlation functions. For n = 4, fif- 
teen components in eq.(4) contribute equally, and thus 
(r?Tf) = (of of) = (T^frJ^) = (4P id - 1)/15. For 
general n, there is a relation, 

(x?xf°) = (<P i;j > - ±) , (5) 

for i ^ j, a ^ 0. 

Although the concept of the high temperature expan- 
sion is quite simple, it needs various devices to calculate 
terms of high orders. We mainly use the finite cluster 
method. 7 ' 8 ^ In this method, series coefficients in the 
thermodynamic limit are obtained exactly by weighted- 
summation of those calculated in finite-size clusters. The 
cluster of size I is defined by 

£-1 

H e :=jJ2 P i,i+i- (6) 

i=l 

To obtain the series expansion of (Pi,j) up to O[(0J) M ], 
we need to calculate Tr[(Ht) m ] and Tr[P itj (H e ) m ] for 
< m < M. In calculating the trace, we adopt a method 
to use combinatorics. 9, 10 ' Then, coefficients are obtained 
as explicit functions of n. Recently a related expansion 
for the specific heat and the uniform susceptibility was 
reported for the Heisenberg model (n = 2). 11 ' We have 
checked that the series obtained by our method repro- 
duces the previous result. 

In using our method, we have a trick to store data 
and reduce computational time and memory. Further, 
we use another trick to reduce the required maximum 
system size £ max . If only the finite cluster method is used, 
the series for (Pj^+i) up to O[(0J) M ] requires £ max — 
M/2 + 1 (M: even) or £ max = (M + l)/2 + 1 (M: odd), 
and the Fourier transform of (Pjj ) requires £ max — M + 
1. However, we reduce £ max by formulating a method 
complementary to the finite cluster method. Namely, we 
can calculate (Pj^+i) up to O[(0J) 22 ] with £ max = 11. 
The Fourier transform of the correlation function (Pij) 
up to O[(0J) 18 ] is calculated with ^ max = 12. As for the 
specific-heat in the n — > oo limit, coefficients of (0J) m 
with odd m are equal to zero, and non-zero coefficients 
are calculated up to O[(0J) 2& }. 



We perform the Pade approximation (PA) and 
the first-order inhomogeneous differential approxima- 
tion 12 ) (FOIDA) to extrapolate the series. In general, the 
FOIDA shows better convergence than the PA. However, 
changing the expansion variable can improve the conver- 
gence of the PA. We change the variable from 0J to 
tanh(c/3J), where c is an adjustable parameter to im- 
prove convergence. 13 -* For small c, the results do not dif- 
fer much from those in 0J. For large c, the convergence 
becomes better. However the tanh variable saturates for 
0J much larger than 1/c, and the approximant tends 
to saturate accordingly. Therefore, we take the smallest 
possible c as long as the convergence remains good. We 
refer to this extrapolation as the PAT hereafter. 

To investigate the specific heat C(T), we first ana- 
lyze the series for the entropy S(T) and then calculate 
C = TdS/dT to use the known fact S(T -> oo) = logn. 
For the PAT, we extrapolate S/T rather than S be- 
cause the exact solution for the ground state 2 ) behaves as 
S oc T at T — > 0. When n = 4, our result for the entropy 
shows a good agreement to that of the quantum Monte 
Carlo method. 6 ^ Although the behavior of the entropy 
seems simple at a sight, the specific heat clearly shows 
a structure composed of a peak at low temperature and 
a shoulder at higher temperature. Figure 1 shows the 
results of the specific heat for various values of n using 
the PAT. At T £ 0.4J, the results of the PAT and the 
FOIDA coincide within the width of the line, and thus 
the estimated error is smaller than the width. The high 
temperature limit is given by C — > (1 — n~ 2 )0 2 J 2 . Note 
that the existence of the shoulder structure is clearly 
seen, which accompanies changes of the sign of the cur- 
vature. In this range the expansion converges well. The 
position of the shoulder seems to be independent of n, 
and the structure remains in the limit n — > oo. We shall 
discuss the origin of this structure later in this paper. 

As n increases, the peak in the specific heat becomes 
sharper and shifts to lower temperature. Such a rapid 
change of a function at very low temperature is difficult 
to extrapolate. Hence the error becomes larger with in- 
creasing n. Particularly, the FOIDA becomes unstable 
at T <, 0.4J; the results of the FOIDA are scattered, 
whereas some results agree with that of the PAT. This is 
because we do not use any information at T — > for the 
FOIDA. As for the PAT, we estimate that the error at 
T ~ 0.2 J is about 5% when n = 6, and smaller than it 
when n < 6. The specific heat at T —> using the PAT 
is by 5% ~ 15% larger than the exact value \n(n — 1)T, 
which is derived from the exact solution. 2 ' 

We remark that the behavior of the specific heat is 
different from that of the l/r 2 -model 14 - ) defined by 

■H 1/r , := jJ2[(NMsm{n(i - j)/N}}- 2 P^. 

i<j 

The inset of Fig.l shows the exact result 15,16 ' for the 
l/r 2 -model, which docs not show the peak-shoulder 
structure in contrast with the nearest-neighbor model. 
Furthermore the T — > limit of the specific heat is given 
by i(n — 1)T, which is 2/n of the nearest-neighbor case. 
We interpret this difference as coming from the differ- 
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ent n-dependence of the spinon velocity. 2 ' 17 ) Namely, the 
spinon velocity in the nearest neighbor model decreases 
as 1/n for large n, while the velocity stays constant in the 
1 / r 2 -model. Such difference is understood intuitively by 
considering the supersymmetric t-J model 18 -* and elimi- 
nating the holes. The spinon velocity is related to the 
Fermi velocity in the non-interacting counterpart and re- 
flects its spectrum. The f/V 2 -modcl has a cusp in the 
bottom of the band, which leads to the finite velocity 
even for vanishing Fermi energy. While the cosine band 
of the nearest-neighbor model gives 1 /n-dependence for 
the Fermi velocity. 
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Fig. 1. Specific heat for n = 2 (dashed line), n = 4 (solid line), 
n = 6 (dash-dotted line) and n = oo (dotted line), in units of the 
gas constant R = 8.31J/(mol • K). The n = oo data are plotted 
only for T > 0.4J. The inset shows the exact solution of the 
inverse-square interaction model. 15,16 ) 



In order to see the temperature-dependent nature of 
the spatial correlation, we show in Fig. 2 the next-nearest- 
neighbor correlation (of of +2 ) for n = 2,4. Here we have 
introduced J' := J/n so as to match the correlation func- 
tion in the high-temperature limit. We compare systems 
with n = 2 

and n = 4 

J' ^2 i T i ■ T 3 + a i ■ °J + ' T j)( (T i ■ • 

The lowest order of the series of (of of +a .) is 0[{(3J) X ] and 
its sign is given by (— l) x . Then we have (of of +2 ) > at 
high temperature. As temperature decreases, (ofof +2 ) 
first grows larger, then starts to decrease and turn to 
negative when n = 4. The reason is the following: Each 
of 15 components tries to align itself antiparallcl. How- 
ever, this cannot be attained simultaneously — e.g. the 
Neel states of both t z and o z makes ferromagnetically 
aligned t z <j z . Therefore, as each antiferromagnetic cor- 
relation becomes larger, each short-range order disturbs 
each other, and thus the correlation turns to a longer 
period at low temperature. 

Such a change of correlation at low temperature 




Fig. 2. Temperature dependence of the next-nearest-neighbor 
correlation (of af^_ 2 ) for n = 2,4. The estimated error is within 
the width of the line. The unit of temperature is given by 
J' = J/n. 



also appears in the Fourier transform of the corre- 
lation function for n > 2. We define 5(g) := 

T, x ( X f X j+x) e ~ lqx (" / ft which docs not depend 
on a or (3 because of the SU(n) symmetry. Thus it is 
reduced to 5(g) = J2 x ( s j s j+x) e ~ lqx when n = 4. A 
naive extrapolation of the series for S(q) shows bad con- 
vergence at several values of q. This is because infor- 
mation at distance x is lacking when cos(xg) ~ 0. To 
avoid this difficulty we extrapolate the series of a com- 
plex function 2 ^ X>() {X^ X^" x )e~ tqx , and take the real 
part of the extrapolated function. This method makes 
use of its imaginary part as well. As a result the conver- 
gence becomes better than the naive extrapolation. Fig- 
ure 3 shows the result at several different temperatures 
for n = 4, 20. The results of the PAT and the FOIDA 
coincide within the width of the line. In the high temper- 
ature limit, only the on-site contribution (XjXj) re- 
mains and S(q) becomes a constant 1/n. As temperature 
decreases, S(tt) increases first, and reaches the maximum 
at T Tmax ~ 0.8 J. With further decrease of temperature, 
S(tt) starts decreasing. Finally the maximum of S(q) 
starts moving from q = it at T mv ~ 0.5 J. Comparing the 
different n results at the same temperature, the peak of 
5(g) around g = ir becomes broader as n increases, es- 
pecially at lower temperature. However, T Tmax and T mv 
seem to be independent of n, even in the n — > oo limit. 

Other numerical studies 5, 6 > for n = 4 show that 5(g) 
at zero temperature has cusps at g = 0,7r/2,7r and has 
the maximum at g = tt/2. For general n, the dominant 
correlation in the ground state 2-4 ) is that of q = 2-k/u, 
and the maximum of 5(g) is expected to reach g = 2n/n 
in the T — > limit. The Fourier transform leading to the 
cusp type behavior requires data of spatial correlations 
for large distance. In Fig. 3 the temperatures are still too 
high to observe this behavior. 

A plausible characteristic temperature of the long- 
range correlation is given by the mean-field approxima- 
tion. Note that the transition temperature in this ap- 
proximation is given by T^f F = 2 J/n, which decreases 
with increasing n. Therefore, if we scale energies by 
T^f F , 5(g) starts to shift toward the ground-state cor- 
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relation at higher temperature as n increases. Namely 
T mv /T^f F increases. 





thcr, interactions making different colors attractive, and 
the same colors repulsive are long-ranged. Then the n- 
site periodic (q = 2n/n) correlation is favorable at low 
temperature. As a result, the same correlation as that 
of the ground state should develop gradually from high 
temperature. The specific heat therefore shows a single 
peak. 

In conclusion, we find by large orders of high- 
temperature expansion a peak-shoulder structure in 
the specific heat which results from different behaviors 
of spatial correlations at high and low temperatures. 
Namely the short-range correlation at high temperature 
has a two-site periodicity, while the ground-state corre- 
lation has the n-site one. We will discuss applications of 
the high-temperature expansion to higher dimensional 
systems in the near future. 
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Fig. 3. Fourier Transform S(q) of the correlation function for (a) 
n = 4 and (b) n = 20 at temperature T = oo (dotted line), 
T = 2J (dash-dotted line), T = 0.8J (dashed line) and T = 0.4J 
(solid line). The estimated error is within the width of the line. 
At T = 0.4J, the maximum is located at q ^ 7r in both (a) and 
(b). 



From the results of the specific heat and the corre- 
lation function, we consider the reason why the struc- 
ture appears in the specific heat. Let us imagine cooling 
down the system gradually from the high temperature 
limit. In the case of the nearest-neighbor interaction, 
short-range correlation with two-site periodicity devel- 
ops first at high temperature. It makes energy gain and 
the specific heat grows as temperature decreases. How- 
ever, since the dominant correlation in the ground state 
is not the q = n component, its energy gain saturates. 
This saturation should cause the shoulder of the spe- 
cific heat. Then, there is another energy gain when the 
system begins to achieve the long-range correlation with 
q = 2ir/n. It should cause the peak at the temperature. 

One may expect that the shift of the peak position of 
the specific heat is mainly due to the shift of T^f F = 2J/n 
to lower temperature with increasing n. However, the 
peak position decreases more rapidly than the 1/n be- 
havior. One of the reasons should be that T^f F is not 
for the q = 2n/n component but for the q = n compo- 
nent. Strong quantum effects are required to realize the 
q = 2ir/n correlation. Accordingly a temperature lower 
than T^f F gives a peak in the specific heat. 

For comparison, let us consider the l/r 2 -model. There 
arc interactions even in the intra-sublattice of the Neel 
state. These make the short-range two-site periodic cor- 
relation less favorable 19 ^ even at high temperature. Fur- 
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